﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using Nature;
using Nature.Numerics;
using System.Diagnostics;
using Nature.Numerics.Basic;

namespace Nature.Toolkit.GEqu
{
    public class UVIdealGasChemicalEquilibriumState : IChemicalEquilibriumState, ISoneFunction
    {
        #region Private Fields
        IChemicalEquilibriumState m_gibbs;
        SoneNewton<UVIdealGasChemicalEquilibriumState> m_newton;
        double m_ro;

        /// <summary>
        /// Mixture internal energy [J/kg]
        /// </summary>
        double m_usum;
        double[] m_x0;
        double[] m_uk;
        double? m_t0, m_p0, m_t, m_p;
        double[] m_tpVector = new double[2];
        bool m_equilibrated = false;
        bool m_compositionInitialized = false;
        #endregion

        #region Private Methods
        void Equilibrate()
        {
            if (m_equilibrated) { return; }
            if (false == m_t0.HasValue) { throw ChemicalEquilibriumStateNotInitializedException.CreateNullTemperatureException(); }
            if (false == m_p0.HasValue) { throw ChemicalEquilibriumStateNotInitializedException.CreateNullPressureException(); }
            if (false == m_compositionInitialized) { throw ChemicalEquilibriumStateNotInitializedException.CreateNullCompositionException(); }

            double wsum = this.Model.CalcMixtureMolarMassFromX(m_x0);

            this.Model.IdealGasSpeciesThermoPotentials.MolarU(m_t0.Value, m_uk);
            m_usum = 0.0;
            for (int isp = 0; isp < Nsp; ++isp) { m_usum += m_x0[isp] * m_uk[isp]; }
            m_usum /= wsum;

            m_ro = wsum * (m_p0.Value / m_t0.Value / Constants.Rgas);

            m_tpVector[0] = m_t ?? m_t0.Value;
            m_tpVector[1] = m_p ?? m_p0.Value;
            object anchor = null;
            m_gibbs.X = m_x0;
            m_newton.Solve(m_tpVector, ref anchor);
            m_t = m_tpVector[0];
            m_p = m_tpVector[1];

            m_equilibrated = true;
        }
        #endregion

        #region Constructors

        public UVIdealGasChemicalEquilibriumState(ModelObject model)
        {
            this.Model = model;
            m_gibbs = new TPIdealGasChemicalEquilibriumState(model);
            m_newton = new SoneNewton<UVIdealGasChemicalEquilibriumState>(this);
            m_newton.SetAbsoluteConvergenceRadius(100.0);
            m_newton.SetRelativeConvergenceRadius(8.0e-1);
            m_uk = new double[model.Nsp];
            m_x0 = new double[model.Nsp];
        }

        #endregion

        #region Public Properties
        public ModelObject Model { get; private set; }

        public int Nsp { get { return this.Model.Nsp; } }
        #endregion

        #region IChEqCalculator Members

        double IChemicalEquilibriumState.T
        {
            get
            {
                Equilibrate();
                return m_t.Value;
            }
            set
            {
                if (m_t0 != value)
                {
                    m_equilibrated = false;
                    m_t0 = value;
                }
            }
        }

        double IChemicalEquilibriumState.P
        {
            get
            {
                Equilibrate();
                return m_p.Value;
            }
            set
            {
                if (m_p0 != value)
                {
                    m_equilibrated = false;
                    m_p0 = value;
                }
            }
        }

        ReadOnlyArray<double> IChemicalEquilibriumState.X
        {
            get
            {
                Equilibrate();
                return m_gibbs.X;
            }
            set
            {
                m_compositionInitialized = true;
                bool changing = false;
                for (int isp = 0; isp < Nsp; ++isp)
                {
                    changing |= (m_x0[isp] != value[isp]);
                    m_x0[isp] = value[isp];
                }
                if (changing)
                {
                    m_equilibrated = false;
                }
            }
        }

        void IChemicalEquilibriumState.DropState()
        {
            m_t0 = m_p0 = null;
            m_t = m_p = null;
            m_compositionInitialized = false;
            m_equilibrated = false;
            m_gibbs.DropState();
        }

        #endregion

        #region ISoneFunction Members

        int ISoneFunction.Nva { get { return 2; } }

        void ISoneFunction.Compute(double[] @x, double[] @y)
        {
            double t = @x[0];
            double p = @x[1];
            m_gibbs.T = t;
            m_gibbs.P = p;
            var _x = m_gibbs.X;
            double wsum = this.Model.CalcMixtureMolarMassFromX(_x);
            this.Model.IdealGasSpeciesThermoPotentials.MolarU(t, m_uk);

            double usum = 0.0;
            for (int isp = 0; isp < Nsp; ++isp) { usum += _x[isp] * m_uk[isp]; }
            usum /= wsum;

            double ro = wsum * (p / t / Constants.Rgas);

            @y[0] = (usum - m_usum);
            @y[1] = (ro - m_ro);

        }

        #endregion
    }
}